

conley_burlig <- function(reg, lat, lon, cutoff, data) {
    
    model.matrix <- model.matrix(reg)
    
    # pull from the model
    n <- nrow(model.matrix)
    k <- ncol(model.matrix) 
    ydata <- model.frame(reg)[,1]
    xdata <- model.matrix
    
    betahat <- solve(t(xdata) %*% xdata) %*% t(xdata) %*% ydata
    e <- ydata - xdata %*% betahat
    
    # grab latitude & longitude
    latdata <-  data[rownames(model.matrix),lat]
    londata <- data[rownames(model.matrix),lon]
    
    # loop over all of the spatial units (aka observations)
    meatWeight <- lapply(1:n, function(i) {
        # turn longitude & latitude into KMs. 1 deg lat = 111 km, 1 deg lon = 111 km* cos(lat)
        lonscale <- cos(latdata[i]*pi / 180) * 111
        latscale <- 111
        # distance --> use pythagorean theorem! who knew that was useful?
        dist <- as.numeric(sqrt((latscale*(latdata[i] - latdata))^2
                                + (lonscale*(londata[i] - londata))^2))
        # set a window var = 1 iff observation j is within cutoff dist of obs i
        window <- as.numeric(dist <= cutoff)
        # this next part is where the magic happens. this thing makes:
        # sum_j(X_iX_j'e_ie_j K(d_{ij})), and we make n of them - one for each i.
        # double transpose here is because R is bad at dealing with 1 x something stuff.
        # we want x_i'; this is an easy way to get it. Now for some dimensions
        # (we want k x k at the end):
        XeeXh <- ((t(t(xdata[i, ])) %*% matrix(1, 1, n) * e[i,]) *
                      # k x 1 1 x n 1 x 1
                      (matrix(1, k, 1) %*% (t(e) * t(window)))) %*% xdata
        # k x 1 1 x n n x k
        return(XeeXh)
    })
    # phew! Now let's make our sandwich. First, the meat = sum_i what we just made
    meat <- (Reduce("+", meatWeight)) / n
    # and the usual bread
    bread <- solve(t(xdata) %*% xdata)
    # mmmm, delicious sandwich
    sandwich <- n* (t(bread) %*% meat %*% bread)
    # se as per usual
    se <- sqrt(diag(sandwich))
    output <- list("betahat" = betahat, 
                   "conleySE" = se)
    return(output)
}
